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Abstract 

We report on a new general relativistic computational model enhancing, in various 
respects, the capability of presently available tools for fitting spectra of X-ray sources. The 
new model is intended for spectral analysis of black-hole accretion discs. Our approach is 
flexible enough to allow easy modifications of intrinsic emissivity profiles. Axial symmetry 
is not assumed, although it can be imposed in order to reduce computational cost of 
data fitting. The main current application of our code is within the XSPEC data-fitting 
package, however, its applicability goes beyond that: the code can be compiled in a stand- 
alone mode, capable of examining time-variable spectral features and doing polarimetry 
of sources in the strong-gravity regime. Basic features of our approach are described in a 
separate paper (Dovciak, Karas & Yaqoob [5]). Here we illustrate some of its applications 
in more detail. We concentrate ourselves on various aspects of line emission and Compton 
' reflection, including the current implementation of the lamp-post model as an example of 

£C) , a more complicated form of intrinsic emissivity. 

m 

O ! 1 Introduction 

Regions of strong gravitational field are most usually explored via X-ray spectroscopy, be- 
cause very hot X-ray emitting material is commonly believed to be present in regions near 
a neutron star surface or a black-hole horizon. Accretion plays crucial role in the process of 
energy liberation and mass accumulation that takes place in this kind of objects 117| . In 
particular, disc-type accretion represents an important mode which is realized under suitable 
circumstances, given by the global geometrical arrangement and local microphysics of the fluid 
medium. A central compact body, undetectable via its own radiation, resides in galactic nuclei 
where it is surrounded by a rather dense population of stars and gaseous environment. Photons 
emerging from the accretion disc and its corona are influenced by gravity of the black hole, so 
they bear various imprints of the gravitational field structure. This concerns especially X-rays 
originating very near the core and, for this reason, spectral analysis with X-ray satellites is 
particularly relevant for astronomical study of strong gravitational fields around black holes. 
For a general discussion, see review articles |1()II27| and further references cited therein. 

In this paper we describe a newly developed computational model aimed for the spectral 
analysis of line profiles and continuum originating in a geometrically thin, planar accretion 
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disc near a rotating (Kerr) black hole. Such analysis has been routinely performed via XSPEC 
package , which performs deconvolution of observed spectra for the effective area and energy 
redistribution of the detector. Previously, several routines were developed and linked with this 
package in order to fit data to a specific model of a black- hole accretion disc ^JED]. However, 
a substantially improved variant of the computational approach has been desirable because 
previous tools have various limitations that may be critical for analyses of present-day and 
forthcoming high-resolution data. 

We describe the layout and usage of the new code and we show some examples and com- 
parisons between new model components. Our present contribution provides information com- 
plementary to the basic description which can be found in ref. [H]- We suggest the reader to 
consult that paper as well as further details in Thesis [5], as they give more examples and 
citations to previous works. Here we concentrate our attention to technical issues of the code 
structure and its performance when computing and fitting spectra. Different perspectives and 
applications of general relativistic computations for black-hole accretion discs have been con- 
sidered by various authors. In particular, it is very useful to consult recent papers of Gierlihski, 
Maciolek-Niedzwiecki & Ebisawa ^31 an d Schnittman & Bcrtschinger |2Ej • Very recently, a new 
independent code has been developed by Beckwith & Done Q]- Their approach also allows to 
study accretion disc spectra including strong gravity effects of a Kerr black hole. This is also 
one of the applications of our code, and so relatively accurate comparisons between both tools 
are possible. We performed several such comparisons and found a very good agreement in the 
shape of predicted line profiles. 

Given a limited space for this contribution, we cannot describe all aspects of the new code: 
capability of the code with respect to timing and polarimetry are discussed elsewhere. Nonethe- 
less, it may be good to bear in mind that such capability has been already implemented and 
tested, taking into account all strong-gravity effects on time-delays and the Stokes parameters. 

2 Transfer functions 

We concentrate on geometrically thin and optically thick accretion discs and we point out that 
general relativity effects can play a role if the configuration is sufficiently dense in a limited 
region, typically a few tens of gravitational radii. Nevertheless, our computational domain 
extends up to about « 10 3 gravitational radii in a non-uniform spatial grid. 

In order to calculate the final spectrum that an observer at infinity measures when local 
emission from the accretion disc is given, one must first specify the intrinsic emissivity in frame 
co-moving with the disc medium and then perform transfer of photons to a distant observer. 
Here we concentrate on the latter part of this task. 

Six functions need to be computed across the source: (i) energy shift affecting the photons 
(i.e. gravitational and Doppler g- factor, needed to account for spectral redistribution), (ii) grav- 
itational lensing (for the evaluation of radiation flux or count rate), (iii) direction of emission 
with respect to the disc normal (for the limb darkening effect), (iv) relative time-delay of the 
light signal (i.e. the mutual delay between photons arriving from different regions of the source, 
needed for the proper account of the light-time effect in timing analysis), (v) change of the 
polarization angle due to photon propagation in the gravitational field (for polarimetry), and 
(vi) azimuthal direction of emission (also for polarimetry). While the first three quantities are 
always necessary, the time-delay factor is required only when local emission is not stationary 
(e.g. the case of orbiting spots) and the change of the polarization angle with the azimuthal 
direction of emission are essential for calculating the overall degree and angle of polarization, 
as observed at infinity. In the adopted approximation of geometrical optics, light rays follow 
null geodesies (in curved space-time) and spectral computations are reduced to a text-book 
problem jS] which, however, may be rather demanding computationally on practical level of 
data fitting. Useful form of light-ray equations and further references can be found in various 
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Figure 1: Six transfer functions are shown in the equatorial plane of a Kerr black hole with 
a = 0.9987 (horizon = 1.05). The black hole rotates counter-clockwise. The observer is 
located upwards at the inclination O = 70°. The values of transfer functions are encoded by a 
colour scale, as indicated above each graph. For mathematical formulae defining these functions 
see eqs. (jA"65f . (|Al^ - (jA"7l"]) and (|A"78)) in Appendix EU 



papers EH]- We summarize the basic equations in Appendix I A. II 

In order for our new model to be fast and practical, we pre-calculated the transfer functions 
for 21 values of the angular momentum of the black hole and 20 values of the inclination angle 
of the observer. The choice of the grid appears sufficiently fine to ensure high accuracy. We 
have stored the transfer functions in the form of tables - i.e. as binary extensions of a FITS 
file. For a technical description of the files layout see Appendix IA.3.11 Values of the transfer 
functions are interpolated when integrating the spectrum for a given angular momentum of the 
black hole and inclination angle of the observer. 

The graphical representation of the tables is shown in Figure ^ Six frames of contour plots 
correspond to individual transfer functions, which are necessary in computations. This figure 
captures equatorial plane for given values of a and 9 . The radius extends up to r — 10 3 r g 
in the tables, but here we show only the central region, r < 6r g , where relativistic effects are 
most prominent. 1 Clock-wise distortion of the contours is due to frame-dragging near a rapidly 
rotating Kerr black hole, and it is clearly visible in the Boyer-Lindquist coordinates here. Notice 

1 Spheroidal coordinates have been employed. We denote r g = GM,/c 2 Ri 1.5 X lO 5 (Af./M0) cm and we use 
geometrized units with c = GM, = 1 hereafter, which means that we scale lengths with M, . Therefore, all 
quantities are dimensionless. 
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that this dramatic distortion appears in the graphical representation only. In order to achieve 
high accuracy of the tables the dragging effect has been largely eliminated in computations by 
means of appropriate transformation of coordinates, as described below. 



3 Photon flux from an accretion disc 

Properties of radiation are described in terms of photon numbers. The source appears as a 
point-like object for a distant observer, so that the observer measures the flux entering the 
solid angle d0 o , which is associated with the detector area dS =D 2 d0 o (see Figure^)- This 
relation defines distance D between the observer and the source. We denote the total photon 
flux received by a detector, 



N*(E) 



where 



dn(E) 
dt dS Q 



dSlN^E/g). 



dm(£i) 
dr dS\ dQ\ 



(1) 



(2) 



is a local photon flux emitted from the surface of the disc, dn(E) is the number of photons with 
energy in the interval {E, E + dE) and g = E/E\ is the redshift factor. The local flux, N\(E\), 
may vary over the disc as well as in time, and it can also depend on the local emission angle. 
This dependency is emphasized explicitly only in the final formula otherwise it is omitted 
for brevity. 






Figure 2: Denomination of various elements of solid angles and areas defined in the text: a) the 
light source appears to the observer to be point-like; b) the light rays received by the detector 
are coming from different parts of the disc (closer view of the disc than in previous figure); c) 
area of a light tube changes as the light rays travel close to the black hole (the disc is edge on). 

The emission arriving at the detector within the solid angle df2 (see Figure^) originates 
from the proper area d<Si on the disc (as measured in the rest frame co-moving with the disc). 
Hence, in our computations we want to integrate the flux contributions over a fine mesh on the 
disc surface. To achieve this aim, we adjust eq. Q to the form 

N Q S (E) = ^JdS^ N x {E/g) 9 2 = ^JdS^±§- NtfS/g) g> . (3) 

Here dSf stands for an element of area perpendicular to light rays corresponding to the solid 
angle df2 at a distance D, dS± is the proper area measured in the local frame of the disc 
and perpendicular to the rays, and dS is the coordinate area for integration. We integrate 
in a two-dimensional slice of a four-dimensional space-time, which is specified by coordinates 
9 = tt/2 and t = t a — At with At being a time delay with which photons from different parts 
of the disc (that lies in the equatorial plane) arrive to the observer (at the same coordinate 
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time t a ). Therefore, let us define the coordinate area by (we employ coordinates t', 0, r, ip with 
H = t - At and At = At(r, 9, p)) 



dS: 



|d 2 S t / 



—d 2 S 6 

dt' btl 



\d 2 S t 



Ad 2 S, 



(4) 



We define the tensor d 2 S a p by two four- vector elements dx^ 1 = (dii,dr, 0, 0) and dx 2 ' = 
(d<2, 0, 0, dip) and by Levi-Civita tensor e a /3 7 s- The time components of these vectors, dti 
and d^2, are such that the vectors dx^ and dx% lie in the tangent space to the above defined 
space-time slice. Then we obtain 



dS 



\g 0d et9afi dx [ ™dx%) | = g 6e \J-\\g^ v \\ dr d ¥> = drdp, 



(5) 



where g^ v is the metric tensor and ||<?^|| is the determinant of the metric. The proper area, 
dS±, perpendicular to the light ray can be expressed covariantly in the following way: 



dS± = 



U a V P d 2 ^ 



a/3 



(6) 



Here, d<Sj_ is the projection of an element of area, defined by d 2 ^^, on a spatial slice of an 
observer with velocity U a and perpendicular to light rays. U a is four- velocity of an observer 
measuring the area dS±, and p@ is four-momentum of the photon. The proper area dS± 
corresponding to the same flux tube is identical for all observers. This means that the last 
equation holds true for any four- velocity U a , and we can express it as 



p> 3 d 2 S a p+p a dS ± -- 
For a = t (note that d 2 5 tr = d 2 S tv = 0) we get 



0. 



dS±_ 
~~dS~ 



1 dS± 



g ee d 2 S, 



w 



Pe_ 

Pt 



t,r,6,<p. 



g 



(7) 



(8) 



In the last equation we used the formula for the cosine of local emission angle /x e , see eq. <|A67|) . 
and the fact that we have chosen such an affine parameter of the light geodesic that p t = — 1. 
From eqs. Q and JHJl we get for the observed flux per unit solid angle 



dn(E) 



N 



dt dn o 

where Nq is a normalization constant and 



dr 



dp Ni(E I g) g I (i e r, 



I 



dSf 
dSl 



(9) 



(10) 



is the lensing factor in the limit D — > oo (the limit is performed while keeping D 2 dfl constant, 
see Figure Et). 

For the line emission, the normalization constant AT is chosen in such a way that the total 
flux from the disc is unity. In the case of a continuum model, the flux is normalized to unity 
at a certain value of the observed energy (typically at E = 1 keV, as in other XSPEC models). 

Finally, the integrated flux per energy bin, AE, is 



AN^(E, AE, t) ■ 



E+AE 



dEN^(E,t) 



dr 



dip 



(E+AE)/g 



E/g 



dEi N x {E h r, p, fi e , t - At) g z r , (II) 
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where At is the relative time delay with which photons arrive to the observer from different parts 
of the disc. The transfer functions g, I, fi e and At are read from the FITS file KBHtablesNN .fits 
described in Appendix IA.3.11 This equation is numerically integrated for a given local flux 
N\(E\, r,ip, He,t — At) in all hereby described new general relativistic non-axisymmetric models. 

Let us assume that the local emission is stationary and the dependence on the axial co- 
ordinate is only through the prescribed dependence on the local emission angle f(fJ> e ) (limb 
darkening/brightening law) together with an arbitrary radial dependence R(r), i.e. 

Ni(E h r, <p, Me, t - At) ee N X (E X ) R(r) /( Me ). (12) 

The observed flux N^(E) is in this case given by 

/oo 
dEkNiiEfiG&Ei), (13) 
-oo 

where 

f'out /■27T 



G(E,E l )=N drR(r) dip f([i c ) g I fj, e r5(E — gE\). (14) 
Jn n Jo 

In this case, the integrated flux can be expressed in the following way: 

i-E+AE 



AN^(E, AE) = / dEN^(E) = 
Je 



E 

E+AE /"Tout r2n foo 

dEN drR(r) d<p f(fi c ) N^E/g) gl Mc r / 6Ey S(E l - E/g) = 

E Jr in JO J -oo 

/•rout i-oo p(E+AE)/E 1 

= N drR(r) dE x N x {E x ) dgF{g), (15) 

Jr ia J-oc Je/Ei 

where we substituted g — E / E\ and 

F(g)= f d V f{^)g 2 l^ c r8{g-g). (16) 
Jo 

Eq. (|15|l is numerically integrated in all axially symmetric models. The function dF(g) = 
dg F(g) has been pre-calculated for several limb darkening/brightening laws /(/i ) and stored 
in separate files, KBHlineNN . f its (see Appendix IA.3.2|I . 

4 Stokes parameters in strong gravity regime 

For polarization studies, Stokes parameters are used. Let us define specific Stokes parameters 
in the following way: 

where I u , Q v , U v and V v are Stokes parameters for light with frequency v, E is energy of 
a photon at this frequency. Further on, we drop the index v but we will always consider 
these quantities for light of a given frequency. We can calculate the integrated specific Stokes 
parameters (per energy bin), i.e. Ai Q , Aq Q , Au and Av Q . These are the quantities that the 
observer determines from the local specific Stokes parameters i\, q\, u\ and v\ on the disc in the 
following way: 

Ai (E,AE) = N J dS j ' dEih(E x ) Fr , (18) 



G 



Aq (E,AE) 
Au Q (E, AE) 
Av (E, AE) 

Here, F = F(r, ip) — g 2 I /i c is & transfer function, iff is the angle by which a vector parallelly 
transported along the light geodesic rotates. We refer to this angle also as a change of the 
polarization angle, because the polarization vector is parallelly transported along light geodesies. 
See Figure |3| for an exact definition of the angle "J. The integration boundaries are the same as 
in eq. (fill) . As can be seen from the definition, the first specific Stokes parameter is equal to 
the photon flux, therefore, eqs. and (|18|l are identical. The local specific Stokes parameters 
may depend on r, ip, fj, e and t — At, which we did not state in the eqs. 

O-GD explicitly for 

simplicity. 

The specific Stokes parameters that the observer measures may vary in time in the case when 
the local parameters also depend on time. In eqs. I|18|) - (|21[l we used a law of transformation of 
the Stokes parameters by the rotation of axes (eqs. (1.185) and (1.186) in 4 ). 

(i) Let three- vectors pa, n>A, n' A and /a be the momen- 
tum of a photon, normal to the disc, projection of the 
normal to the plane perpendicular to the momentum 
and a vector which is parallelly transported along the 
geodesic (as four- vector) , respectively; 

(ii) let iff a be an angle between n' A and /a; 

(iii) let the quantities in (i) and (ii) be evaluated at the 
disc for A—\ with respect to the local rest frame co- 
moving with the disc, and at infinity for A = 2 with 
respect to the stationary observer at the same light 
geodesic; 

(iv) then the change of polarization angle is defined as 
* = * 2 - 



Figure 3: Definition of the change of polarization angle iff. 

An alternative way for expressing polarization of light is by using the degree of polarization 
P and two polarization angles Xo and £ D , defined by 

= y/ql+ul + vl/i , (22) 
tan2xo = u /q , (23) 
sin 2^ = Vo/y/ql+ul + v*. (24) 

5 New model for XSPEC 

We have developed several general relativistic models for line emission and Compton reflection 
continuum. The line models are supposed to be more accurate and versatile than the LAOR 
model JHj) and substantially faster than the kerrspec model [201 - Several models of intrinsic 
emissivity were employed, including the lamp-post model |22| . Among other features, these 
models allow various parameters to be fitted such as the black-hole angular momentum, observer 
inclination, accretion disc size and some of the parameters characterizing disc emissivity and 



= N j dS J dE\[qx{E x )cos2iff -u\{E x )sw.2^]Fr, (19) 
= N J dS J dEi [qi(Ei) sin 2* + Ui(-Ei) cos 2*] Fr , (20) 
= N J dS J dEi vi(Ei) Fr . (21) 
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primary illumination properties. They also allow a change in the grid resolution and, hence, to 
control accuracy and computational speed. Furthermore, we developed very general convolution 
models. All these models are based on pre-calculated tables described in SectionEland thus the 
geodesies do not need to be calculated each time one integrates the disc emission. These tables 
are calculated for the vacuum Kerr space-time and for a Keplerian co-rotating disc plus matter 
that is freely falling below the marginally stable orbit. The falling matter has the energy and 
angular momentum of the matter at the marginally stable orbit. It is possible to use different 
pre-calculated tables if they are stored in a specific FITS file (see Appendix lA.3.1l for its detailed 
description) . 

There are two types of new models. The first type integrates the local disc emission in 
both of the polar coordinates on the disc and thus enables one to choose non-axisymmetric 
area of integration. This option is useful for example when computing spectra of spots or 
partially obscured discs. One can also choose the resolution of integration and thus control 
the precision and speed of the computation. The second type of models is axisymmetric - the 
axially dependent part of the emission from rings is pre-calculated and stored in a FITS file 
(the function dF(g) = dgF(g) from lltjfl is integrated for different radii with the angular grid 
having 20 000 points). These models have less parameters that can be fitted and thus are less 
flexible even though more suited to the standard analysis approach. On the other hand they 
are fast because the emission is integrated only in one dimension (in the radial coordinate of 
the disc). It may be worth emphasizing that the assumption about axial symmetry concerns 
only the form of intrinsic cmissivity of the disc, which cannot depend on the polar angle in this 
case, not the shape of individual light rays, which is always complicated near a rotating black 
hole. 
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unit 
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maximum value 
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GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM/c 2 


400. 


0. 


999. 


zshif t 




0. 


-0.999. 


10. 


ntable 




0. 


0. 


99. 



Table 1: Common parameters for all models. 

There are several parameters and switches that are common for all new models (see Table^): 

a/M - specific angular momentum of the Kerr black hole in units of GM/c (M is the mass of 

the central black hole), 
theta_o - the inclination of the observer in degrees, 

rin-rh - inner radius of the disc relative to the black-hole horizon in units of GM/c 2 , 
ms - switch for the marginally stable orbit, 

rout-rh - outer radius of the disc relative to the black- hole horizon in units of GM/c 2 , 
zshif t - overall redshift of the object, 

ntable - number of the FITS file with pre-calculated tables to be used. 

The inner and outer radii are given relative to the black-hole horizon and, therefore, their 
minimum value is zero. This becomes handy when one fits the a/M parameter, because the 
horizon of the black hole as well as the marginally stable orbit changes with a/M, and so the 
lower limit for inner and outer disc edges cannot be set to constant values. The ms switch 
determines whether we intend to integrate also emission below the marginally stable orbit. If 
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its value is set to zero and the inner radius of the disc is below this orbit then the emission 
below the marginally stable orbit is taken into account, otherwise it is not. 

The ntable switch determines which of the pre-calculated tables should be used for intrinsic 
emissivity. In particular, ntable = for KBHtablesOO . f its (KBHlineOO . f its), ntable = 1 
for KBHtablesOl . f its (KBHlineOl . f its), etc., corresponding to non-axisymmetric (axisym- 
metric) models. 



parameter 


unit 


default value 


minimum value 


maximum value 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


Stokes 




0. 


0. 


6. 



Table 2: Common parameters for non-axisymmetric models. 

The following set of parameters is relevant only for non-axisymmetric models (see Table |2J): 

phi - position angle of the axial sector of the disc in degrees, 
dphi - inner angle of the axial sector of the disc in degrees, 
nrad - radial resolution of the grid, 

division - switch for spacing of radial grid (0 - equidistant, 1 - exponential), 
nphi — axial resolution of the grid, 

smooth - switch for performing simple smoothing (0 - no, 1 - yes), 
Stokes - switch for computing polarization (see Table I3J). 

The phi and dphi parameters determine the axial sector of the disc from which emission comes 
(see Figure 0J. The nrad and nphi parameters determine the grid for numerical integration. If 
the division switch is zero, the radial grid is equidistant; if it is equal to unity then the radial 
grid is exponential (i.e. more points closer to the black hole). 

If the smooth switch is set to unity then a simple smoothing is applied to the final spectrum. 
Here A^(£,) = [jV n (£i_i) + 2N?(E } ) + N^(E i+1 )}/4. 

If the Stokes switch is different from zero, then the model also calculates polarization. Its 
value determines, which of the Stokes parameters should be computed by xspec, i.e. what will 
be stored in the output array for the photon flux photar; see Table 01 (If Stokes ^ then 
a new ascii data file stokes.dat is created in the current directory, where values of energy 
E together with all Stokes parameters i, q, u, v, P, x[deg] and £[deg] are stored, each in one 
column.) 

A realistic model of polarization has been currently implemented only in the KYLlCR model 
(see Section below). In other models, a simple assumption is made - the local emission is 
assumed to be linearly polarized in the direction perpendicular to the disc (i.e. q\ = i\ = N\ and 
u\ = Vi = 0). In all models (including KYLlCR) there is always no final circular polarization 
(i.e. v = £ = 0), which follows from the fact that the fourth local Stokes parameter is zero in 
each model. 

6 Models for a relativistic spectral line 

Three general relativistic line models are included in the new set of xspec routines - non- 
axisymmetric Gaussian line model kygIline, axisymmetric Gaussian line model kygline and 
fluorescent lamp-post line model kyfIll. 
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Figure 4: Segment of a disc from which emission comes (view from above). 



value 


photon flux array photar contains^ 





i = 


I/E, where I is the first Stokes parameter (intensity) 


1 


q = 


- Q/E, where Q is the second Stokes parameter 


2 


u = 


- U/E, where U is the third Stokes parameter 


3 


V = 


= V/E, where V is the fourth Stokes parameter 


4 


de£ 


;ree of polarization, P = q 2 + u 2 + v 2 /i 


5 


anj 


;le x[deg] of polarization, tan 2% = u/q 


6 


anj 


de C[dcg] of polarization, sin 2^ = vj \J q 2 + u 2 + v 2 



1 the photar array contains values described in the table and multiplied by 
width of the corresponding energy bin 
t E is energy of observed photons 



6.1 Non-axisymmetric Gaussian line model KYGlLINE 

The KYGlLINE model computes the integrated flux from the disc according to eq. 1)1 ljl. It 
assumes that the local emission from the disc is 



The local emission is assumed to be a Gaussian line with its peak flux depending on the radius 
as a broken power law. The line is defined by nine points equally spaced with the central point 
at its maximum. The normalization constant No in 1)11)1 is such that the total integrated flux 
of the line is unity. The parameters defining the Gaussian line are (see Table HJ|: 

Erest - rest energy of the line in keV, 

sigma - width of the line in eV, 

alpha - radial power-law index for the outer region, 

beta - radial power-law index for the inner region, 

rb - parameter defining the border between regions with different power-law indices, 



Table 3: Definition of the Stokes parameter. 




for r > rb 



(25) 




for r <r h . (26) 
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jump - ratio between flux in the inner and outer regions at the border radius, 
limb - switch for different limb darkening/brightening laws. 

There are two regions with different power-law dependences with indices alpha and beta. The 
power law changes at the border radius where the local emissivity does not need to be 
continuous (for jump 7^ 1). The rb parameter defines this radius in the following way: 

r h = rb x r ms for rb > , (27) 
r b = -rb + r h for rb < , (28) 

where r ms is the radius of the marginally stable orbit and rh is the radius of the horizon of the 
black hole. 



parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 


- 


1. 


0. 


1. 


rout-rh 


GM/c 2 


400. 


0. 


999. 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 


*Erest 


keV 


6.4 


1. 


99. 


* sigma 


eV 


2. 


0.01 


1000. 


*alpha 




3. 


-20. 


20. 


*beta 




4. 


-20. 


20. 


*rb 




0. 


0. 


160. 


*jump 




1. 


0. 


Ie6 


*limb 




-1. 


-10. 


10. 


Stokes 




0. 


0. 


6. 



Table 4: Parameters of the non-axisymmetric Gaussian line model KYGlLlNE. Model parame- 
ters that are not common for all non-axisymmetric models are denoted by asterisk. 

The function /(^t e ) = /(cos 6 C ) in Ij25(l and (|26|l describes the limb darkening/brightening 
law, i.e. the dependence of the local emission on the local emission angle. Several limb darken- 
ing/brightening laws are implemented: 

/Owe) = 1 for limb = , (29) 

/Oite) = 1 + 2.06^ for limb = -1 , (30) 

/(/x e ) = lna + A^ 1 ) for limb = -2, (31) 

/(/i e ) = Mc imb for limb ^0,-1, -2. (32) 

Eq. I)29|l corresponds to the isotropic local emission, eq. (|30[) corresponds to limb darkening 
in an optically thick electron scattering atmosphere (used by Laor |18l 1191 and eq. Q31J1 
corresponds to limb brightening predicted by some models of a fluorescent line emitted by an 
accretion disc due to X-ray illumination |131 ITS] . 
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There is also a similar model kyg2line present among the new xspec models, which is 
useful when fitting two general relativistic lines simultaneously. The parameters are the same 
as in the KYGlLlNE model except that there are two sets of those parameters describing the 
local Gaussian line emission. There is one more parameter present, ratio21, which is the ratio 
of the maximum of the second local line to the maximum of the first local line. Polarization 
computations are not included in this model. 

6.2 Axisymmetric Gaussian line model KYGLINE 

This model uses eq. I|15|) for computing the disc emission with local flux being 

Ni(Ei) = 5(Ei - Erest) , (33) 
R(r) = r ~ al P ha . (34) 

The function dF(g) = dgF(g) in l|16fl was pre-calculated for three different limb darken- 
ing/brightening laws l)29|) - (|31|l and stored in corresponding FITS files KBHlineOO . f its - 
KBHline02 . f its. The local emission is a delta function with its maximum depending on the 
radius as a power law with index alpha and also depending on the local emission angle. The 
normalization constant Nq in (|15|l is such that the total integrated flux of the line is unity. 



parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM/c 2 


400. 


0. 


999. 


zshif t 




0. 


-0.999 


10. 


ntable 




1. 


0. 


99. 


*Erest 


keV 


6.4 


1. 


99. 


*alpha 




3. 


-20. 


20. 



Table 5: Parameters of the axisymmetric Gaussian line model kygline. Model parameters 
that are not common for all axisymmetric models are denoted by asterisk. 

There are less parameters defining the line in this model than in the previous one (see 
Table EJ): 

Erest - rest energy of the line in keV, 
alpha - radial power-law index. 

Note that the limb darkening/brightening law can be chosen by means of the ntable switch. 

This model is much faster than the non-axisymmetric kygIline model. Although it is not 
possible to change the resolution grid on the disc, it is hardly needed because the resolution 
is set to be very large, corresponding to nrad = 500, division = 1 and nphi = 20 000 in 
the kygIline model, which is more than sufficient in most cases. (These values apply if the 
maximum range of radii is selected, i.e. rin=0, ms=0 and rout=999; in case of a smaller range 
the number of points decreases accordingly.) This means that the resolution of the kygIline 
model is much higher than what can be achieved with the laor model, and the performance is 
still very good. 

6.3 Non-axisymmetric fluorescent lamp-post line model KYF1 LL 

The line in this model is induced by the illumination of the disc from the primary power-law 
source located on the axis at height above the black hole. This model computes the final 
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spectrum according to eq. (|llf> with the local photon flux 



PhoIndex-l sin #Ld6> L 
9h 



• dr 



x exp 



1000 



E\ — Ere st 



2height 
height 2 + c 

2" 



/(Mi) Me 



V2 



sigma 



(35) 



Here, <7l is ratio of the energy of a photon received by the accretion disc to the energy of the 
same photon when emitted from a source on the axis, $l is an angle under which the photon 
is emitted from the source (measured in the local frame of the source) and /Zi = cos Si is the 
cosine of the incident angle (measured in the local frame of the disc) - see Figure 03 All of 
these functions depend on height above the black hole at which the source is located and on 
the rotational parameter a/M of the black hole. Values of <?l, 6*l and /ii for a given height and 
rotation are read from the lamp-post tables lamp, fits (see Appendix I A . 3 . 3JI . At present, only 
tables for a/M = 0.9987492 (i.e. for the horizon of the black hole r h = 1.05) and height = 
2,3, 4, 5, 6, 8, 10, 12, 15, 20, 30, 50, 75 and 100 are included in lamp. fits, therefore, the a/M 
parameter is used only for the negative values of height (see below). 

The factor in front of the function /(/Ltj,/i G ) gives the radial dependence of the disc emis- 
sivity, which is different from the assumed broken power law in the kygIline model. For the 
derivation of this factor, which characterizes the illumination from a primary source on the axis 
see Appendix I A. 21 



parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM/c 2 


400. 


0. 


999. 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 


*PhoIndex 




2. 


1.5 


3. 


*height 


GM/c 2 


3. 


-20. 


100. 


*Erest 


keV 


6.4 


1. 


99. 


*sigma 


eV 


2. 


0.01 


1000. 


Stokes 




0. 


0. 


6. 



Table 6: Parameters of the fluorescent lamp-post line model kyfIll. Model parameters that 
are not common for all non-axisymmetric models are denoted by asterisk. 

The function /(/x;, /i e ) is a coefficient of reflection. It depends on the incident and reflection 
angles. Although the normalization of this function also depends on the photon index of the 
power-law emission from a primary source, we do not need to take this into account because 
the final spectrum is always normalized to unity. Values of this function are read from a pre- 
calculated table which is stored in fluorescent J.ine .fits file (see an d Appendix I A . 3 . 4p . 
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The local emission 1)350 is denned in nine points of local energy E\ that are equally spaced 
with the central point at its maximum. The normalization constant No in the formula is 
such that the total integrated flux of the line is unity. The parameters defining local emission 
in this model are (see Table EJ: 

Pholndex - photon index of primary power-law illumination, 

height - height above the black hole where the primary source is located for height > 0, and 

radial power-law index for height < 0, 
Erest - rest energy of the line in keV, 
sigma - width of the line in eV. 

If positive, the height parameter works as a switch - the exact value present in the tables 
lamp, fits must be chosen. If the height parameter is negative, then this model assumes that 
the local emission is the same as in the KYGIline model with the parameters alpha = —height, 
rb — and limb = —2 (Pholndex parameter is unused in this case). 




Figure 5: Reflection models: a) lamp-post model; b) diffuse corona model. 



7 Compton reflection models 

We have developed two new relativistic continuum models - lamp-post Compton reflection 
model kylIcr and the kyhIrefl model which is a relativistically blurred hrefl model that 
is already present in XSPEC. Both of these models are non-axisymmetric. 

7.1 Non-axisymmetric lamp-post Compton reflection model KYLlCR 

The emission in this model is induced by the illumination of the disc from the primary power- 
law source located on the axis at height above the black hole. As in every non-axisymmetric 
model the observed spectrum is computed according to eq. (|ll|l . The local emission is 

N X {E{) = r he± z ht f(E i; fi c ) for height < . (37) 

For the definition of ^l, #l and /ii see Section f6.3l and Appendix IA.3.31 where pre-calculated 
tables of these functions in lamp, fits are described. 

The function f(E\; fii, /J. e ) gives dependence of the locally emitted spectrum on the angle of 
incidence and the angle of emission, assuming a power-law illumination. This function depends 
on the photon index Pholndex of the power-law emission from a primary source. Values of 
this function for various photon indices of primary emission are read from the pre-calculated 
tables stored in ref spectra. fits (see Appendix IA.3.5|I . These tables were calculated by 
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the Monte Carlo simulations of Compton scattering |2J. At present, tables for Pho Index = 
1.5, 1.6, . . . , 2.9, 3.0 and for local energies in the range from 2 keV to 300 keV are available. 
The normalization constant No in eq. (|ll(l is such that the final photon flux at an energy of 
3 keV is equal to unity, which is different from what is usual for continuum models in XSPEC 
(where the photon flux is unity at 1 keV). The choice adopted is due to the fact that current 
tables in ref spectra. fits do not extend below 2 keV. 

The function f(E x ; /i c ), which is used for negative height, is an averaged function f(E x ; /x e ) 
over jU; 

f(E x ;» e ) = [ d/ii/GEu/ii,/^). (38) 



The local emission (|37ll can be interpreted as emission induced by illumination from clouds 
localized near above the disc rather than from a primary source on the axis (see Figure EJ). In 
this case photons strike the disc from all directions. 

For positive values of height the KYLlCR model includes a physical model of polarization 
based on Rayleigh scattering in single scattering approximation. The specific local Stokes 
parameters describing local polarization of light are 

i x (E x ) = h T +I * N X (E X ), (39) 



(h 


+ Ir) 


h 


-I, 


(h 


+ Ir) 




u 



91 = rhrr\ N ^ El ^ ( 4 °) 

{h + It) 

ux{E x ) = j^-jjN x {E x ), (41) 
v x {E x ) = 0, (42) 

where functions I x , I r and U determine the angular dependence of the Stokes parameters in the 
following way 



Mc 2 (l + /if) + 2(1 - _ tf) _ ^^(l-^Xl-^) COS ($ c - $i) 

-/i2(l- M ?)cos2($ e -$i), (43) 
I t = l + M f + (i_ M 2) C0S 2($ c -$ i ), (44) 



h 



U = -4^(1-^X1-^) sin($ c - $ i )-2// c (l-//| ! )sin2($ -$i), (45) 

Here 4> e and $i are the azimuthal emission and the incident angles in the local rest frame 
co- moving with the accretion disc (see Appendixes IA.1I and IA.2I for their definition). For the 
derivation of these formulae see the definitions (1.147) and eqs. (X.172) in 4. We have omitted a 
common multiplication factor, which would be cancelled anyway in eqs. (|39J) - H42|) . The symbol 
( ) in definitions of the local Stokes parameters means value averaged over the difference of 
the azimuthal angles <i> — $i. We divide the parameters by (I\ + I t ) because the function 
f(E x ; fi x ,fi c ), and thus also the local photon flux N\(E\), is averaged over the difference of the 
azimuthal angles. 

The parameters defining local emission in this model are (see Tabled: 
Pholndex - photon index of primary power-law illumination, 

height - height above the black hole where the primary source is located for height > 0, and 

radial power-law index for height < 0, 
line - switch whether to include the iron lines (0 - no, 1 - yes), 
E_cut - exponential cut-off energy of the primary source in keV. 
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parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM /c 


400. 


0. 


999. 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 


*PhoIndex 




2. 


1.5 


3. 


*height 


GM/c 1 


3. 


-20. 


100. 


^1 -i no 

x j. lie 




n 

u. 


n 


1 

1 . 


*E_cut 


keV 


300. 


1. 


1000. 


Stokes 




0. 


0. 


6. 



Table 7: Parameters of the lamp-post Compton reflection model kylIcr. Model parameters 
that are not common for all non-axisymmetric models are denoted by asterisk. 

Tables ref spectra. fits for the function /(i?i; /Zj, /Lte) also contain the emission in the iron 
lines Ka and K/3. The two lines can be excluded from computations if the line switch is set 
to zero. The E_cut parameter sets the upper boundary in energies where the emission from a 
primary source ceases to follow a power-law dependence. If the E_cut parameter is lower than 
both the maximum energy of the considered dataset and the maximum energy in the tables for 
f(E\; fii, fie) in ref spectra. fits (300 keV), then this model is not valid. 

7.2 Non-axisymmetric Compton reflection model KYHlREFL 

This model is based on an existing multiplicative hrefl model in combination with the pow- 
ERLAW model, both of which are present in XSPEC. Local emission in <|11[1 is the same as 
the spectrum given by the model hrefl*powerlaw with the parameters thetamin = and 
thetamax = 90 with a broken power-law radial dependence added: 

Ni(Ei) = r~ alplla hrefl*powerlaw for r > r h , (46) 
iVi(-Ei) = jumpr£ eta ~ alpha r~ beta HREFL*POWERLAW for r < r h . (47) 

For a definition of the boundary radius by the rb parameter see eqs. (|27|) (I28|) . and for a 
detailed description of the hrefl model see |Hj and the XSPEC manual. The kyhTrefl model 
can be interpreted as a Compton-reflection model for which the source of primary irradiation 
is near above the disc, in contrast to the lamp-post scheme with the source on the axis (see 
Figure EJ). The approximations for Compton reflection used in hrefl (and therefore also in 
kyhIrefl) are valid below ~ 15 keV in the disc rest-frame. The normalization of the final 
spectrum in this model is the same as in other continuum models in XSPEC, i.e. photon flux is 
unity at the energy of 1 keV. 

The parameters defining the local emission in kyhIrefl (see Table [SJ arc 

Pho Index - photon index of the primary power-law illumination, 
alpha - radial power-law index for the outer region, 
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parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM/c z 


400. 


0. 


999. 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 


*PhoIndex 




1. 


0. 


10. 


*alpha 




3. 


-20. 


20. 


*beta 




4. 


on 
-20. 


OA 

20. 


*rb 




0. 


0. 


160. 


*jump 




1. 


0. 


le6 


*Feabun 




1. 


0. 


200. 


*FeKedge 


keV 


7.11 


7.0 


10. 


Jib CI IdL 




i 

i . 


n 

u. 


1UUU. 


*covf ac 




1. 


0. 


1000. 


Stokes 




0. 


0. 


6. 



Table 8: Parameters of the reflection kyhIrefl model. Model parameters that are not common 
for all non-axisymmetric models are denoted by asterisk. 

beta - radial power-law index for the inner region, 

rb - parameter defining the border between regions with different power-law indices, 
jump - ratio between flux in the inner and outer regions at the border radius, 
Feabun - iron abundance relative to solar, 
FeKedge - iron K-edge energy, 

Escf rac - fraction of the direct flux from the power-law primary source seen by the observer, 
covf ac - normalization of the reflected continuum. 

8 General relativistic convolution models 

We have also produced two convolution-type models, kyIconv and kyconv, which can be 
applied to any existing XSPEC model for the intrinsic X-ray emission from a disc around a 
Kerr black hole. We must stress that these models are substantially more powerful than the 
usual convolution models in xspec (these are commonly defined in terms of one-dimensional 
integration over energy bins). Despite the fact that our convolution models still use the stan- 
dard xspec syntax in evaluating the observed spectrum (e.g. kyconv(powerlaw)), our code 
accomplishes a more complex operation. It still performs ray-tracing across the disc surface so 
that the intrinsic model contributions are integrated from different radii and azimuths on the 
disc. 

There are several restrictions that arise from the fact that we use existing xspec models: 
- by local xspec models only the energy dependence of the photon flux can be defined, 
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- only a certain type of radial dependence of the local photon flux can be imposed - we 
have chosen to use a broken power-law radial dependence, 

- there is no azimuthal dependence of the local photon flux, except through limb darkening 



- local flux depends on the binning of the data because it is defined in the centre of each 
bin, a large number of bins is needed for highly varying local flux. 

For emissivities that cannot be defined by existing xspec models, or where the limitations 
mentioned above are too restrictive, one has to add a new user-defined model to xspec (by 
adding a new subroutine to xspec). This method is more flexible and faster than convolution 
models (especially when compared with the non-axisymmetric one), and hence it is recom- 
mended even for cases when these prefabricated models could be used. In any new model for 
xspec one can use the common ray-tracing driver for relativistic smearing of the local emission: 
ide for non-axisymmetric models and idre for axisymmetric ones. For a detailed description 
see Appendixes lA.4.1l and lA.4.2l 

8.1 Non-axisymmetric convolution model KYClONV 

The local emission in this model is computed according to eq. (| 1 1 1> with the local emissivity 
equal to 



For a definition of the boundary radius by the rb parameter see eqs. ()27|) - (I28|) and for 
definition of different limb darkening laws /(/x e ) see eqs. H29|) - ()32|) . The local emission is given 
by the model in the centre of energy bins used in XSPEC with the broken power-law radial 
dependence and limb darkening law added. Apart from the parameters of the model, the local 
emission is defined also by the following parameters (see Table EJ: 

normal - switch for the normalization of the final spectrum, 
= - total flux is unity (used usually for the line), 

> - flux is unity at the energy = normal keV (used usually for the continuum), 

< - flux is not normalized, 
ne_loc - number of points in the energy grid where the local photon flux is defined, 
alpha - radial power-law index for the outer region, 
beta - radial power-law index for the inner region, 

rb - parameter defining the border between regions with different power-law indices, 
jump - ratio between the flux in the inner and outer regions at the border radius, 
limb - switch for different limb darkening/brightening laws. 

The local emission in each ky model has to be defined either on equidistant or exponential (i.e. 
equidistant in logarithmic scale) energy grid. Because the energy grid used in the convolution 
model depends on the binning of the data, which may be arbitrary, the flux has to be rebinned. 
It is always rebinned into an exponentially spaced energy grid in ky convolution models. The 
ne_loc parameter defines the number of points in which the rebinned flux will be defined. 

8.2 Axisymmetric convolution model KYCONV 

The local emission in this model is computed according to eq. 115JI with the local emissivity 
equal to 



law. 




f(/J.e) MODEL for r < r h . 



(48) 
(49) 



Ni(Ei) = MODEL, 

R(r) = r~ al P ha . 



(50) 
(51) 
(52) 
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parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 




1. 


0. 


1. 


rout-rh 


GM /c 


400. 


0. 


999. 


phi 


deg 


0. 


-180. 


180. 


dphi 


deg 


360. 


0. 


360. 


nrad 




200. 


1. 


10000. 


division 




1. 


0. 


1. 


nphi 




180. 


1. 


20000. 


smooth 




1. 


0. 


1. 


*normal 




1. 


-1. 


100. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 


*ne_loc 




100. 


3. 


5000. 


* alpha 




3. 


-20. 


20. 


*beta 




4. 


on 

-20. 


OA 

20. 


*rb 


'ms 


0. 


0. 


160. 


*jump 




1 

1 . 


n 
u. 


reo 


*limb 




0. 


-10. 


10. 


Stokes 




0. 


0. 


6. 



Table 9: Parameters of the non-axisymmetric convolution model KYClONV. Model parameters 
that are not common for all non-axisymmetric models are denoted by asterisk. 

Except for the parameters of the model, the local emission is defined also by the following 
parameters (see Table HHjl : 

alpha - radial power-law index, 

ne_loc - number of points in energy grid where the local photon flux is defined, 
normal - switch for the normalization of the final spectrum, 
= - total flux is unity (used usually for the line), 

> - flux is unity at the energy = normal keV (used usually for the continuum), 
< - flux is not normalized. 

Note that the limb darkening/brightening law can be chosen through the ntable switch. This 
model is much faster than the non-axisymmetric convolution model KYClONV. 

9 Examples and comparisons 

In our new models we have concentrated ourselves mainly on two components that contribute 
to the X-ray spectra of active galactic nuclei and X-ray binaries with black-hole candidates 
- spectral line emission and its relativistic broadening, and the Compton reflection from an 
illuminated disc. Two basic types of illumination have been considered - the disc illuminated 
either from every direction by a nearby diffuse corona above the disc, or from a particular direc- 
tion by a small source placed on the axis above the black hole (see Figure [SJ. The illumination 
in the former case decreases with radius as a power law. Hence, this model is characterized 
by the radial power-law index a. On the other hand, the illumination anisotropy in the latter 
(lamp-post) model depends mainly on the position of a primary source of emission characterized 
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parameter 


unit 


default value 


minimum value 


maximum value 


a/M 


GM/c 


0.9982 


0. 


1. 


theta_o 


deg 


30. 


0. 


89. 


rin-rh 


GM/c 2 


0. 


0. 


999. 


ms 


~ 


1. 


0. 


1. 


rout-rh 


GM/c 2 


400. 


0. 


999. 


zshif t 




0. 


-0.999 


10. 


ntable 




0. 


0. 


99. 






Q 
O. 


-ZU. 


zu. 


*ne_loc 




100. 


3. 


5000. 


^normal 




1. 


-1. 


100. 



Table 10: Parameters of the axisymmetric convolution model kyconv. Model parameters that 
are not common for all axisymmetric models are denoted by asterisk. 




0.1 1 10 0.1 1 10 

Energy [keV] Energy [keV] 

Figure 6: An example of a line profile originating from a disc in equatorial plane of a Kerr black 
hole (a = 0.9987 GM/c 2 , i.e. — 1.05 GM/c 2 ) due to the illumination from a primary source 
on the axis. The kyfIll model was used. Left: Dependence of the line profile on the height (in 
GM/c 2 ) of a primary source with photon index r = 2. Right: Dependence of the line profile 
on the photon index of the primary emission with a source at height 3 GM/c 2 above the black 
hole. 

by height h where it is located. In both cases it is assumed that the primary emission has a 
continuum power-law shape which can be characterized by a photon index T (Pholndex). 

The emission from the disc depends on quite a number of parameters. It is influenced by the 
mass M and the rotation a of the central black hole, by the area from which the emission from 
disc comes (defined by inner radius r- m , outer radius r out and azimuthal segment with boundaries 
at ip and ip + Aip), by the inclination O of the observer, by the radial power-law index a, by 
the photon index T and by the position h of a primary source. Limb darkening/brightening 
law (dependence on the local emission angle) is another important factor that determines the 
final spectrum we observe. 

Here, we will show several examples of emission for the lamp-post fluorescent line model and 
for the reflection models. For other examples and comparisons, see the accompanying paper 
(Dovciak, Karas & Yaqoob 9 ). In all figures in this section, we assumed the inclination angle 
o = 30°, the rotational parameter a = 0.9987 GM/c 2 , and an emitting ring extending from 
r in = r ms to r out = AQQGM/c 2 . 

In Figure we demonstrate that the broad iron emission lines due to illumination from the 
source placed on the axis depend heavily on the height where the "lamp" is located (left), as 
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Figure 7: Comparison of the general relativistic kyhIrefl model with non-relativistic 
hrefl(powerlaw). The relativistic blurring of the iron edge is clearly visible. The power-law 
index of the primary source is PhoIndex=2 (left) and PhoIndex=2.6 (right). 
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Figure 8: General relativistic lamp-post Compton reflection model kylIcr with (dashed) and 
without (solid) iron lines Ka and K/3. The emission from the disc is induced by illumination 
from a primary source placed 2GM/c 2 (left) and 100 GM/c 2 (right) above the black hole. 



well as on the photon index of the primary emission (right). It can be seen that the intrinsic 
width of the line (2 eV in this example) is much less than its subsequent relativistic broadening, 
and the local profile (assumed to be Gaussian) is thus smeared in the final spectrum. These 
graphs correspond to the iron Ka line with the rest energy of 6.4 keV. 

Relativistic effects are demonstrated also in Figure [7] where the non-relativistic reflection 
model hrefl(powerlaw) is compared with our relativistic kyhIrefl. Blurring of the iron 
edge is clearly visible. Here, we set the radial power-law index a = 1 in kyhIrefl. Other 
parameters defining these models were set to their default values. 

Examples of the Compton reflection emission component of the spectra with and without 
the fluorescent Ka and K/3 lines are shown in Figure |H1 It can be seen that originally narrow 
lines can contribute substantially to the continuum component. 

We compare the two new relativistic reflection models kylIcr and kyhIrefl in Figure 
Note that the kylIcr model is valid only above approximately 2 keV and the kyhIrefl model 
only below approximately 15keV. 
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Figure 9: Comparison of the two new general relativistic Compton reflection models kylIcr 
and kyhIrefl. The lamp-post kylIcr model is characterized by the height h above the disc 
where a primary source of emission is placed, the reflection kyhIrefl model is characterized by 
the radial power-law index a. Left: h = 2 GM/c 2 , a = 3.4 . Right: h = 100 GM/c 2 , a = 1.5 . 

10 Conclusions 

In this paper we described the main features of the newly developed set of routines. We 
have concentrated ourselves on various technical issues connected with fitting X-ray spectra 
using our model. In particular, we described several variants of the code which are suited for 
modelling relativistic spectral components originating in a Keplerian disc near a rotating black 
hole. Both axially symmetric and non-axisymmetric models were discussed. For further details 
and for exemplary analysis of XMM-Newton satellite data we refer to Dovciak et al. 0| ■ 

Our package offers a number of applications which could not be examined in the limited 
space of the present paper. In particular, timing analysis can be performed with the code, 
but we defer detailed description of this capability to subsequent papers. Also, we have only 
briefly touched the possibility of polarimetric analysis, which offers great possibilities for future 
X-ray spectroscopy but goes beyond routine capabilities of devices installed onboard present 
day satellites. Additional emissivity laws can be easily adopted. This can be achieved either by 
using the convolution component or by adding a new user-defined model. The latter method is 
more flexible and faster, and hence recommended. In both approaches, the ray-tracing routine 
is linked and used for relativistic blurring. 

As general motivation for developing this project further, we remind the reader that various 
disc-like structures are almost ubiquitous in objects where the fluid orbits around and inflows 
onto a compact body. The central mass, M,, can vary by many orders of magnitude in different 
objects, and its value provides the basic classification for black-hole sources. Physical character- 
istics of accretion discs also scale roughly with M, . Indeed, accretion discs around supermassive 
black holes in active galactic nuclei and quasars share some properties with circumstellar discs 
in close binary systems, e.g. cataclysmic variable stars and microquasars. However, there are 
important distinctions between the two kinds of objects which prohibit any simple scaling (for 
example, galactic nuclear discs tend to be colder and less dense compared to circumstellar discs). 
In both cases there is strong evidence suggesting that some spectral components (namely, the 
iron Ka line emission) originate, at least in part, within ~ 10 gravitational radii of a central 
black-hole. 

The central compact body governs gravitational field in which the medium of an accretion 
flow evolves. Since we consider a general relativistic description of the gravitational field, the 
rotation of the central body should not be ignored. An angular momentum is actually one of 
the model parameters which could in principle be measured by means of spectral analysis of 
observed radiation. 
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A Appendix 

A.l Summary of equations 

Before writing equations for the transfer functions let us summarize basic formulae defining the 
Kerr space-time, light geodesies and disc's motion. We remind the reader that units GM, = 
c = 1 are used (M, is the mass of the central black hole). 
The Kerr metric in Boyer-Lindquist coordinates is 
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where p 2 = r 2 + a 2 cos 2 #, A = r 2 — 2r + a 2 and A = (r 2 + a 2 ) 2 — Aa 2 sin 2 #. We assume < a < 1 
everywhere in this paper. 

The four-momentum of the photons emitted from the disc is (see e.g. and |23| ) 

pi = [a{l-a) + (r 2 +a 2 )(r 2 +a 2 ~al)/A]/r 2 , (A53) 

p r c = R, gn {{r 2 +a 2 -al) 2 -A[{l-a) 2 +q 2 ]} 1 ' 2 /r 2 , (A54) 

pi = -q/r 2 , (A55) 

pt = [l-a + a{r 2 + a 2 -al)/A]/r 2 . (A56) 

Here I = asin# and q 2 — fi 2 + cos 2 (a 2 — a 2 ) are constants of motion with a and (3 being 
impact parameters measured perpendicular and parallel, respectively, to the spin axis of the 
black hole projected onto the observer's sky. Here we define a to be positive when photon 
travels in the direction of the four- vector ^ at infinity, and [3 to be positive if it travels in the 

direction of — ^ at infinity. Furthermore, we have denoted sign of the radial component of the 
momentum by R sgn . We have chosen an affine parameter along light geodesies in such a way 
that the conserved energy is normalized to — p c t = — Pot = 1- 

The Keplerian velocity of the co-rotating disc above the marginally stable orbit is 

U< = /' 2 + aV¥ , (A57) 
T v r 2 — 3r + 2ay / r 

U r = 0, (A58) 

U 6 = 0, (A59) 

U v = 1 (A60) 

yjr(r 2 - 3r + 2a^F) 

We assume that the matter in the disc below the marginally stable orbit conserves its specific 
energy and its specific angular momentum, i.e. Ut{r < r ms ) = Ut{r ms ) and U v {r < r ms ) = 
U v (r ms ). We get the radial component U r (r < r ms ) from the normalization of the four- velocity, 
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■1, and from the fact that the disc rotates in the equatorial plane even below the 



marginally stable orbit, i.e. U (r < r ms ) = 0. 

In our calculations we use the following local orthonormal tetrad connected with the matter 
in the disc 



e (t)M 

e (0)M 



fg^U r U t 



Vi + u r u, 

(0, 0, V5^, 0) 
~AU V 



, V9rr(l + U r U r ), 0, 



y/l + U r U r 



VI + U r U r 



, 0, 0, 



A U 



(A61) 
(A62) 
(A63) 

(A64) 



The gravitational and Doppler shift (g-factor) is defined as the ratio of the energy of a 
photon received by an observer at infinity to the local energy of the same photon when emitted 
from an accretion disc 

v Q Pot 1 

.9 = — 



(A65) 



Here v Q and v c denote frequency of the observed and emitted photons, respectively. 

We define lensing as the ratio of the area at infinity perpendicular to the light rays through 
which photons come to the proper area on the disc perpendicular to the light rays and corre- 
sponding to the same flux tube 



dSf 
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The four- vectors Y ei and Y e2 are transported along the geodesic according to the equation of the 
geodesic deviation from infinity where they are unit, space-like and perpendicular to each other 
and to the four-momentum of light. In (|A66|I we have denoted the magnitude of a four- vector 
by 1 1 || and the scalar product of two four- vectors by < , > . 
The cosine of the local emission angle is 



/i c = cos <5 C = 



p ca n" 

~p ef iUH 



(A67) 



where n a — —efg) is the normal to the disc with respect to the observer co-moving with the 
matter in the disc. 

The relative time delay At is the Boyer-Lindquist time which elapses between the emission 
of a photon from the disc and its reception by an observer (plus a certain constant so that the 
delay is finite close to the black hole) . We have integrated the equation of the geodesies in Kerr 
ingoing coordinates and thus we have calculated the delay in the Kerr ingoing time coordinate 
Aix- The Boyer-Lindquist time coordinate can be obtained from the Kerr ingoing one by the 
following equation 

2r 



dt = dt 
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with r± = 1 ± Vl — a? . Then we define the delay as 
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(A68) 

for a < 1 , (A69) 
(A70) 



There is a minus in front of the brackets because the direction of integration is from infinity 
(represented by r Q = 10 11 in our computations) to the disc. 
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The change of the polarization angle is (see UJ, 0) 



where 



Y 

tan* = — , (A71) 



X = -(a - asin0 o )rci - pKz , (A72) 
Y = (a — asm0 o )K2 — 0K\ , (A73) 



(A74) 



with Ki and K2 being components of the complex constant of motion K pw (see |2(j| ) 



Kl = arpif -rlapl-ir 2 + a 2 )pt]f -r(r 2 +a 2 )p e c r , (A75) 
K2 = -rplf + rlpl-ap^r + arpir. (A76) 

Here the polarization vector is a four-vector corresponding to the three-vector f% from 
Figure which is chosen in such a way that it is a unit vector parallel with n' x (i.e. \&i = 0) 

We define the azimuthal emission angle as the angle between the projection of the three- 
momentum of the emitted photon into the disc (in the local rest frame co-moving with the disc) 
and the radial tetrad vector: 

$ c = -R° gn arccos (^=|) + \ ■ (A78) 

where Rg gn is positive if the emitted photon travels outwards (pi r%> > 0) and negative if it travels 

inwards (j>[^ < 0) in the local rest frame of the disc. 

We conclude this section by the relationship between the Boyer-Lindquist coordinate tp and 
the Kerr ingoing coordinate ifK, which we use when we interpolate between the pre-calculated 
tables 

CL T — V 1 

<P = <PK H In for a < 1 , (A79) 

r + — r_ r — r_ 

ip = ip K — for a = 1 . (A80) 

r — 1 

A. 2 Local emission in lamp-post models 

The local emission from a disc is proportional to the incident illumination from a power-law 
primary source placed on the axis at height h above the black hole. To calculate the incident 
illumination we need to integrate the geodesies from the source to the disc. 

The four-momentum of the incident photons which were emitted by a primary source and 
which are striking the disc at radius r is (see e.g. and |23| ) 

p\ = l + 2/r + 4/A, (A81) 

Pi = K & A{r 2 +a 2 ) 2 -Ha 2 + q l)] 1/2 /r\ (A82) 

Pi = Qh/r 2 , (A83) 

pf = 2o/(rA), (A84) 
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where g£ = sin 2 #L {h 2 + a 2 ) 2 /Al — a 2 is Carter's constant of motion with Al = h 2 — 2h + a 2 , 
and with the angle of emission 6>l being the local angle under which the photon is emitted 
from a primary source (it is measured in the rest frame of the source). We define this angle 

by tan6>L = —Pl/Plj where p[} = e Lu ano - Pl^ ~ PL e Lu with p£ and e-j® being the 
four-momentum of emitted photons and the local tetrad connected with a primary source, 
respectively. The angle is 0° when the photon is emitted downwards and 180° if the photon is 
emitted upwards. 

We denoted the sign of the radial component of the momentum by Rg gn . We have chosen 
such an affine parameter for the light geodesic that the conserved energy of the light is — p\ t — 
—Pht — 1- The conserved angular momentum of incident photons is zero (Zl = 0). 

The gravitational and Doppler shift of the photons striking the disc which were emitted by 
a primary source is 

■9l = — = = TTT- ■ ( A85 ) 



Here v\ and denote the frequency of the incident and emitted photons, respectively and U£ 
is a four-velocity of the primary source with the only non-zero component [/£ = yj h2 ^ 2 ■ 
The cosine of the local incident angle is 

Mi = |cos^| = ^^, (A86) 

where n a = — is normal to the disc with respect to the observer co-moving with the matter 
in the disc. 

We further define the azimuthal incident angle as the angle between the projection of the 
three-momentum of the incident photon into the disc (in the local rest frame co-moving with 
the disc) and the radial tetrad vector, 



— 1 Pia e ( v ) 



$i = -Rl gn arccos ( n = — ] + - , (A87) 



where Rg gn is positive if the incident photon travels outwards (pj > 0) and negative if it travels 

inwards (p^ < 0) in the local rest frame of the disc. 

In lamp-post models the emission of the disc will be proportional to the incident radiation 
Nf(Ei) which comes from a primary source 

N i s (E l ) = N?(E L )^. (A88) 

Here N[^(El) = A^l r is an isotropic and stationary power-law emission from a primary 
source which is emitted into a solid angle d^L and which illuminates local area d<Si on the disc. 
The energy of the photon striking the disc (measured in the local frame co-moving with the 
disc) will be redshifted 

E l= g L E L . (A89) 



The ratio df^/dSi is 
where (see eqs. (0 and J7J) 



d» L = d» L dS = sinfl L dfl L dy> dS 

dS*i dS dSi drdip dSi ' [ ) 



dS = drd<p=\d 2 S t \ = ~g ee ^dS ± = ^dS±. (A91) 

Pi PI 
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Here we used the same space-time slice as in the discussion above eq. (@J and thus the element 
d 2 S a p is defined as before, see eq. Note that here the area dS± is defined by the incident 
flux tube as opposed to dS± in eq. (jHJl where it was defined by the emitted flux tube. The 
coordinate area dS corresponds to the proper area dS± which is perpendicular to the incident 
light ray (in the local rest frame co- moving with the disc). The corresponding proper area 
(measured in the same local frame) lying in the equatorial plane will be 

dS l = \d 2 S (t P | = |ef t) e< 9 > " d 2 V | = \g gg 1/2 U» d 2 S, e \ = 

= -9ee ,2P -^^=g^ 2 % 9lj dS^. (A92) 

Pi Pi 

Here we have used eq. (@J for the tetrad components of the element d 2 S a p, eqs. (© and l|A85(l . 

It follows from eqs. I)A88|1 - (|A92|) that the incident radiation will be again a power law with 
the same photon index T as in primary emission 

N i s (E l ) = N oi E 1 - T , (A93) 

with the normalization factor 



The emission of the disc due to illumination will be proportional to this factor. 

A.3 Description of FITS files 

A. 3.1 Transfer functions in KBHtablesNN . f its 

The transfer functions are stored in the file KBHtablesNN . fits as binary extensions and 
parametrized by the value of the observer inclination angle 9 a and the horizon of the black 
hole rh. We found parametrization by fh more convenient than using the rotational parameter 
a, although the latter choice may be more common. Each extension provides values of a par- 
ticular transfer function for different radii, which are given in terms of r — r^, and for the Kerr 
ingoing axial coordinates ifK ■ Values of the horizon , inclination O , radius r — and angle 
tfK, at which the functions are evaluated, are defined as vectors at the beginning of the FITS 
file. 



The definition of the file KBHtablesNN . fits: 

0. All of the extensions defined below are binary. 

1. The first extension contains six integers defining which of the functions is present in the 
tables. The integers correspond to the delay, g-factor, cosine of the local emission angle, 
lensing, change of the polarization angle and azimuthal emission angle, respectively. Value 
means that the function is not present in the tables, value 1 means it is. 

2. The second extension contains a vector of the horizon values in GM/c 2 (1.00 < n, < 2.00). 

3. The third extension contains a vector of the values of the observer's inclination angle 8 Q 
in degrees (0° < 9 a < 90°, 0° - axis, 90° - equatorial plane). 

4. The fourth extension contains a vector of the values of the radius relative to the horizon 
r — rh in GMj c 2 . 

5. The fifth extension contains a vector of the values of the azimuthal angle tpK in radians 
(0 < <^k < 27r). Note that <pk is a Kerr ingoing axial coordinate, not the Boyer-Lindquist 
one! 

6. All the previous vectors have to have values sorted in an increasing order. 
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7. In the following extensions the transfer functions are defined, each extension is for a par- 
ticular value of rh and O . The values of and 9 Q are changing with each extension in 
the following order: 

rh[i] x o [i], 

r h [l] x o [2], 
r h [l] x O [3], 



r h [2] x 9 [1], 
r h [2] x o [2], 
r h [2] x O [3], 

8. Each of these extensions has the same number of columns (up to six). In each column, 
a particular transfer function is stored - the delay, ^-factor, cosine of the local emission 
angle, lensing, change of the polarization angle and azimuthal emission angle, respectively. 
The order of the functions is important but some of the functions may be missing as defined 
in the first extension (see 1. above). The functions are: 

delay - the Boyer-Lindquist time in GM/c 3 that elapses between the emission of a photon 
from the disc and absorption of the photon by the observer's eye at infinity plus a 
constant, 

g-factor - the ratio of the energy of a photon received by the observer at infinity to the 
local energy of the same photon when emitted from an accretion disc, 

cosine of the emission angle - the cosine of the local emission angle between the emitted 
light ray and local disc normal, 

lensing - the ratio of the area at infinity perpendicular to the light rays through which 
photons come to the proper area on the disc perpendicular to the light rays and 
corresponding to the same flux tube, 

change of the polarization angle in radians - if the light emitted from the disc is linearly 
polarized then the direction of polarization will be changed by this angle at infinity - 
counter-clockwise if positive, clockwise if negative (we are looking towards the coming 
emitted beam); on the disc we measure the angle of polarization with respect to the 
"up" direction perpendicular to the disc with respect to the local rest frame; at 
infinity we also measure the angle of polarization with respect to the "up" direction 
perpendicular to the disc - the change of the polarization angle is the difference 
between these two angles, 

azimuthal emission angle in radians - the angle between the projection of the three-mo- 
mentum of an emitted photon into the disc (in the local rest frame co-moving with 
the disc) and the radial tetrad vector. 

For mathematical formulae defining the functions see eqs. (|A65I . (|A66|I - (|A67() . (|A69|I - 
(|ATTl and l|A"78)l in Appendix |A~T1 

9. Each row corresponds to a particular value of r — rj, (see 4. above). 

10. Each element corresponding to a particular column and row is a vector. Each element of 
this vector corresponds to a particular value of ifK (see 5. above). 

We have pre-calculated three sets of tables - KBHtablesOO.f its, KBHtables50 . f its and 
KBHtables99 . f its. All of these tables were computed for an accretion disc near a Kerr black 
hole with no disc corona present. Therefore, ray-tracing in the vacuum Kerr space-time could 
be used for calculating the transfer functions. When computing the transfer functions, it was 
supposed that the matter in the disc rotates on stable circular (free) orbits above the marginally 
stable orbit. The matter below this orbit is freely falling and has the same energy and angular 
momentum as the matter which is on the marginally stable orbit. 
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The observer is placed in the direction ip = ir/2. The black hole rotates counter-clockwise. 
All six functions are present in these tables. 

Tables are calculated for these values of the black-hole horizon: 
KBHtablesOO.fits: 1.00, 1.05, 1.10, 1.15, 1.90, 1.95, 2.00 (21 elements), 
KBHtables50.fits: 1.00, 1.10, 1.20, 1.90, 2.00 (11 elements), 
KBHtables99.f its: 1.05 (1 element), 
and for these values of the observer's inclination: 

- KBHtablesOO . fits: 0.1, 1, 5, 10, 15, 20, . . . , 80, 85, 89 (20 elements), 

- KBHtables50.f its: 0.1, 1, 10, 20, 80, 89 (11 elements), 
KBHtables99.fits: 0.1, 1, 5, 10, 15, 20, 80, 85, 89 (20 elements). 

The radii and azimuths at which the functions are evaluated are same for all three tables: 

- radii r — are exponentially increasing from to 999 (150 elements), 

- values of the azimuthal angle <^k are equidistantly spread from to 27r radians with a much 
denser cover "behind" the black hole, i.e. near tp K = 1.5-7T (because some of the functions are 
changing heavily in this area for higher inclination angles, o > 70°) (200 elements). 

A. 3. 2 Tables in KBHlineNN . fits 

Pre-calculated functions dF(g) = dg F(g) defined in eq. (fTB)l are stored in FITS files KBHlineNN .fits. 
These functions are used by all axisymmetric models. They are stored as binary extensions and 
they are parametrized by the value of the observer inclination angle O and the horizon of the 
black hole r^. Each extension provides values for different radii, which are given in terms of 
r — rh, and for different g-factors. Values of the g-factor, radius r — 7*h, horizon rh, and incli- 
nation 6 0l at which the functions are evaluated, are defined as vectors at the beginning of the 
FITS file. 

The definition of the file KBHlineNN. fits: 

0. All of the extensions defined below are binary. 

1. The first extension contains one row with three columns that define bins in the ^-factor: 

- integer in the first column defines the width of the bins (0 - constant, 1 - exponentially 
growing), 

- real number in the second column defines the lower boundary of the first bin (minimum 
of the g-factor), 

- real number in the third column defines the upper boundary of the last bin (maximum 
of the g-factor). 

2. The second extension a contains vector of the values of the radius relative to the horizon 
r - r h in GM/c 2 . 

3. The third extension contains a vector of the horizon values in GM/c 2 (1.00 < rh < 2.00). 

4. The fourth extension contains a vector of the values of the observer's inclination angle O 
in degrees (0° < o < 90°, 0° - axis, 90° - equatorial plane). 

5. All the previous vectors have to have values sorted in an increasing order. 

6. In the following extensions the functions dF(g) are defined, each extension is for a par- 
ticular value of rh and O . The values of and O are changing with each extension in 
the same order as in tables in the KBHtablesNN . f its file (see the previous section, point 
7.). Each extension has one column. 

7. Each row corresponds to a particular value of r — (see 2. above). 

8. Each element corresponding to a particular column and row is a vector. Each element 
of this vector corresponds to a value of the function for a particular bin in the ^-factor. 
This bin can be calculated from number of elements of the vector and data from the first 
extension (see 1. above). 
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We have pre-calculated several sets of tables for different limb darkening/brightening laws 
and with different resolutions. All of them were calculated from tables in the KBHtablesOO .fits 
(see the previous section for details) and therefore these tables are calculated for the same values 
of the black-hole horizon and observer's inclination. All of these tables have equidistant bins 
in the g-factor which fall in the interval (0.001, 1.7). Several sets of tables are available: 

KBHlineOO . f its for isotropic emission, see eq. (|29[1 . 

KBHlineOl . f its for Laor's limb darkening, see eq. (|30|l . 

- KBHline02 . f its for Haardt's limb brightening, see eq. (51) ■ 

All of these tables have 300 bins in the ^-factor and 500 values of the radius r — which are 
exponentially increasing from to 999. We have produced also tables with a lower resolution - 
KBHline50.f its, KBHline51 . f its and KBHline52 . f its with 200 bins in the g-factor and 300 
values of the radius. 

A. 3. 3 Lamp-post tables in lamp, fits 

This file contains pre-calculated values of the functions needed for the lamp-post model. It is 
supposed that a primary source of emission is placed on the axis at a height h above the Kerr 
black hole. The matter in the disc rotates on stable circular (free) orbits above the marginally 
stable orbit and it is freely falling below this orbit where it has the same energy and angular 
momentum as the matter which is on the marginally stable orbit. It is assumed that the corona 
between the source and the disc is optically thin, therefore ray-tracing in the vacuum Kerr 
space-time could be used for computing the functions. 

There are five functions stored in the lamp. fits file as binary extensions. They are 
parametrized by the value of the horizon of the black hole r^, and height h, which are de- 
fined as vectors at the beginning of the FITS file. Currently only tables for = 1.05 (i.e. 
a = 0.9987492) and h = 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 30, 50, 75 and 100 are available. The 
functions included are: 

- angle of emission in degrees - the angle under which a photon is emitted from a primary 

source placed at a height h on the axis above the black hole measured by a local stationary 
observer (0° - a photon is emitted downwards, 180° - a photon is emitted upwards), 

- radius - the radius in GM/c 2 at which a photon strikes the disc, 

- g-factor - the ratio of the energy of a photon hitting the disc to the energy of the same photon 

when emitted from a primary source, 

- cosine of the incident angle - an absolute value of the cosine of the local incident angle be- 

tween the incident light ray and local disc normal, 

- azimuthal incident angle in radians - the angle between the projection of the three-momentum 

of the incident photon into the disc (in the local rest frame co-moving with the disc) and 
the radial tetrad vector. 
For mathematical formulae defining the functions see eqs. l|A85|l - ljA87(l in Appendix IA.2I 
The definition of the file lamp, fits: 

0. All of the extensions defined below are binary. 

1. The first extension contains a vector of the horizon values in GM/c 2 , though currently 
only FITS files with tables for one value of the black-hole horizon are accepted (1.00 < 
r h < 2.00). 

2. The second extension contains a vector of the values of heights ft of a primary source in 
GM/c 2 . 

3. In the following extensions the functions are defined. Each extension is for a particular 
value of rh and h. The values of rh and h are changing with each extension in the following 



30 



order: 
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n,[2] 
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X 



X 



X 



Ml], 

ft[2], 
fc[3], 



4. Each of these extensions has five columns. In each column, a particular function is stored 
- the angle of emission, radius, ^-factor, cosine of the local incident angle and azimuthal 
incident angle, respectively. Extensions may have a different number of rows. 

A. 3. 4 Coefficient of reflection in f luorescent_line . f its 

Values of the coefficient of reflection /(/ii,/i c ) for a fluorescent line are stored for different 
incident and reflection angles in this file. For details on the model of scattering used for 
computations see |21| . It is assumed that the incident radiation is a power law with the photon 
index T — 1.7. The coefficient does not change its angular dependences for other photon indices, 
only its normalization changes (see Figure 14 in 22)- The FITS file consists of three binary 
extensions: 

- the first extension contains absolute values of the cosine of incident angles, 

- the second extension contains values of the cosine of reflection angles, 

- the third extension contains one column with vector elements, here values of the coefficient 
of reflection are stored for different incident angles (rows) and for different reflection angles 
(elements of a vector). 

A. 3. 5 Tables in ref spectra . fits 

The function f(E\; fi c ) which gives dependence of a locally emitted spectrum on the angle of 
incidence and angle of emission is stored in this FITS file. The emission is induced by a power- 
law incident radiation. Values of this function were computed by the Monte Carlo simulations 
of Compton scattering, for details see [2] • The reflected radiation depends on the photon index 
r of the incident radiation. There are several binary extensions in this FITS file: 

- the first extension contains energy values in keV where the function f(E\; ^;,/i e ) is com- 
puted, currently the interval from 2 to 300 keV is covered, 

- the second extension contains the absolute values of the cosine of the incident angles, 

- the third extension contains the values of the cosine of the emission angles, 

- the fourth extension contains the values of the photon indices T of the incident power law, 
currently tables for T — 1.5, 1.6, . . . , 2.9 and 3.0 are computed, 

in the following extensions the function f(E\\ /d, fj, c ) is defined, each extension is for a 
particular value of F; here values of the function are stored as a vector for different 
incident angles (rows) and for different angles of emission (columns), each element of this 
vector corresponds to a value of the function for a certain value of energy. 

A. 4 Description of the integration routines 

Here we describe the technical details about the integration routines, which act as a common 
driver performing the ray-tracing for various models of the local emission. The description of 
non-axisymmetric and axisymmetric versions are both provided. An appropriate choice depends 
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on the form of intrinsic emissivity. Obviously, non-axisymmetric tasks are computationally more 
demanding. 

A. 4.1 Non-axisymmetric integration routine ide 

This subroutine integrates the local emission and local Stokes parameters for (partially) po- 
larized emission of the accretion disc near a rotating (Kerr) black hole (characterized by the 
angular momentum a) for an observer with an inclination angle 9 . The subroutine has to be 
called with ten parameters: 

ide (ear , ne , nt , f ar , qar , uar , var , ide_par am , emiss ivity , ne JLoc) 

ear - real array of energy bins (same as ear for local models in xspec), 

ne - integer, number of energy bins (same as ne for local models in xspec), 

nt - integer, number of grid points in time (nt = 1 means stationary model), 

far(ne,nt) - real array of photon flux per bin (same as photar for local models in XSPEC but 

with the time resolution), 
qar(ne,nt) - real array of the Stokes parameter Q divided by the energy, 
uar (ne ,nt) - real array of the Stokes parameter U divided by the energy, 
var (ne ,nt) - real array of the Stokes parameter V divided by the energy, 
ide_param - twenty more parameters needed for the integration (explained below), 
emissivity - name of the external emissivity subroutine, where the local emission of the disc 

is defined (explained in detail below), 
ne_loc - number of points (in energies) where local photon flux (per keV) in the emissivity 

subroutine is defined. 

The description of the ide_param parameters follows: 
ide_param(l) - a/M - the black-hole angular momentum (0 < a/M < 1), 

ide_param(2) theta_o - the observer inclination in degrees (0° - pole, 90° - equatorial 
plane) , 

ide_param(3) - rin-rh - the inner edge of the non-zero disc emissivity relative to the black- 
hole horizon (in GM/c 2 ), 

ide_param (4) - ms - determines whether we also integrate emission below the marginally stable 
orbit; if its value is set to zero and the inner radius of the disc is below the marginally 
stable orbit then the emission below this orbit is taken into account, if set to unity it is 
not, 

ide_param(5) - rout-rh - the outer edge of the non-zero disc emissivity relative to the black- 
hole horizon (in GM/c 2 ), 

ide_param(6) - phi - the position angle of the axial sector of the disc with non-zero emissivity 
in degrees, 

ide_param(7) - dphi - the inner angle of the axial sector of the disc with non-zero emissivity 

in degrees (dphi < 360°), 
ide_param(8) - nrad - the radial resolution of the grid, 

ide_param(9) - division - the switch for the spacing of the radial grid (0 - equidistant, 1 - 
exponential) , 

ide_param(10) - nphi - the axial resolution of the grid, 

ide_param(ll) - smooth - the switch for performing simple smoothing (0 - no, 1 - yes), 
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ide_param(12) - normal - the switch for normalizing the final spectrum, 
if = - total flux is unity (usually used for the line), 

if > - flux is unity at the energy = normal keV (usually used for the continuum), 
if < - final spectrum is not normalized, 
ide_param(13) - zshift - the overall redshift of the object, 

ide_param(14) - ntable - tables to be used, it defines a double digit number NN in the name 
of the FITS file KBHtablesNN.fits containing the tables (0 < ntable < 99), 

ide_param(15) - edivision - the switch for spacing the grid in local energies (0 - equidistant, 
1 - exponential), 

ide_param(16) - periodic - if set to unity then local emissivity is periodic if set to zero it is 

not (need not be set if nt = 1), 
ide_param(17) - dt - the time step (need not be set if nt = 1), 

ide_param(18) polar - whether the change of the polarization angle and/or azimuthal emis- 
sion angle will be read from FITS tables (0 - no, 1 - yes), 
ide_param(19) - rO-rh and 

ide_param(20) - phiO - in dynamical computations the initial time will be set to the time 
when photons emitted from the point [rO, phiO] on the disc (in the Boyer-Lindquist 
coordinates) reach the observer. 

The ide subroutine needs an external emissivity subroutine in which the local emission and 
local Stokes parameters arc defined. This subroutine has twelve parameters: 

emissi vity (ear J.oc, ne_loc,nt, far _Loc,qar_loc,uar_Loc,var_loc,r, phi, cosine, 
phiphoton , f irst_emis ) 

ear_loc(0:ne_Loc) - real array of the local energies where local photon flux far_loc is de- 
fined, with special meaning of ear_loc(0) - if its value is larger than zero then the local 
emissivity consists of two energy regions where the flux is non-zero; the flux between these 
regions is zero and ear_loc(0) defines the number of points in local energies with the 
zero local flux, 

ne_loc - integer, the number of points (in energies) where the local photon flux (per keV) is 
defined, 

nt - integer, the number of grid points in time (nt = 1 means stationary model), 
far_loc(0:ne_Loc,nt) - real array of the local photon flux (per keV) - if the local emissivity 

consists of two separate non-zero regions (i.e. ear_loc(0) > 0) then f ar_loc (0 , it) is the 

index of the last point of the first non-zero local energy region, 
qar_loc (ne_loc ,nt) - real array of the local Stokes parameter Q divided by the local energy, 
uar_loc(ne_loc,nt) - real array of the local Stokes parameter U divided by the local energy, 
var_loc(ne_loc,nt) - real array of the local Stokes parameter V divided by the local energy, 
r - the radius in GM/c 2 where the local photon flux f ar_loc at the local energies ear_loc is 

demanded 

phi - the azimuth (the Boyer-Lindquist coordinate (p) where the local photon flux f ar_loc at 
the local energies ear_loc is demanded, 

cosine - the cosine of the local angle between the emitted ray and local disc normal, 

phiphoton - the angle between the emitted ray projected onto the plane of the disc (in the 
local frame of the moving disc) and the radial component of the local tetrad (in radians), 

f irst_emis - boolean, TRUE if we enter the emissivity subroutine from the ide subroutine for 
the first time, FALSE if this subroutine was already evaluated during the present run. This 
distinction is convenient to initialize some variables when calling the emissivity subroutine 
for the first time (e.g. calculation of the falling spot trajectory can be performed in this 
place) . 
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A. 4. 2 Axisymmetric integration routine idre 

This subroutine integrates the local axisymmetric emission of an accretion disc near a rotating 
(Kerr) black hole (characterized by the angular momentum a) for an observer with an inclination 
angle O . The subroutine has to be called with eight parameters: 

idre (ear ,ne , phot ar , idre_param, cmodel ,ne _loc , ear_Loc , f ar_loc) 

ear - real array of energy bins (same as ear for local models in xspec), 

ne - integer, the number of energy bins (same as ne for local models in xspec), 

photar - real array of the photon flux per bin (same as photar for local models in xspec), 

idre_param - ten more parameters needed for the integration (explained below), 

cmodel - 32-byte string with a base name of a FITS file with tables for axisymmetric emission 

(e.g. "KBHline" for KBHlineNN . fits), 
ne_loc - the number of points (in energies) where the local photon flux (per keV) is defined in 

the emissivity subroutine, 
ear_loc - array of the local energies where the local photon flux f ar_loc is defined 
f ar_loc - array of the local photon flux (per keV). 

The description of the idre_param parameters follows: 

idre_param(l) - a/M - the black-hole angular momentum (0 < a/M < 1), 
idre_param(2) - theta_o the observer inclination in degrees (0° - pole, 90° - equatorial 
plane) , 

idre_param(3) - rin-rh - the inner edge of the non-zero disc emissivity relative to the black- 
hole horizon (in GM/c 2 ), 

idre_param(4) - ms - determines whether we also integrate emission below the marginally 
stable orbit; if its value is set to zero and the inner radius of the disc is below the 
marginally stable orbit then the emission below this orbit is taken into account, if set to 
unity it is not, 

idre_param(5) - rout-rh - the outer edge of the non-zero disc emissivity relative to the black- 
hole horizon (in GM/c 2 ), 
idre_param(6) - smooth - the switch for performing simple smoothing (0 - no, 1 - yes), 
idre_param(7) - normal - the switch for normalizing the final spectrum, 
if = - total flux is unity (usually used for the line), 

if > - flux is unity at the energy = normal keV (usually used for the continuum), 
if < - final spectrum is not normalized, 
idre_param(8) - zshift - the overall redshift of the object, 

idre_param(9) - ntable - tables to be used, it defines a double digit number NN in the name 

of the FITS file (e.g. in KBHlineNN . fits) containing the tables (0 < ntable < 99), 
idre_param(10) - alpha - the radial power-law index. 

The idre subroutine does not need for its operation any external emissivity subroutine. 
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